Thermodynamic Scaling of Diffusion in Supercooled Lennard-Jones Liquids 
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The manner in which the intermolecular potential u{r) governs structural relaxation in liquids is 
a long standing problem in condensed matter physics. Herein we show that diffusion coefficients 
for simulated Lennard-Jones m-6 liquids (8 < m < 36) in normal and moderately supercooled 
states are a unique function of the variable p'' /T, where p is density and T is temperature. The 
scaling exponent 7 is a material specific constant whose magnitude is related to the steepness of the 
repulsive part of u{r), evaluated around the distance of closest approach between particles probed 
in the supercooled regime. Approximations of u{r) in terms of inverse power laws are also discussed. 

PACS numbers: 61.43.Fs, 61.20.Lc, 64.70.Pf, 61.20.Ja 



Establishing a quantitative connection between the re- 
laxation properties of a liquid and the interactions among 
its constituent molecules is the sine qua non for funda- 
mental understanding and prediction of the dynamical 
properties. The supercooled regime is of particular in- 
terest, since both intermolecular forces and steric con- 
straints (excluded volume) exert significant effects on the 
dynamics. This makes temperature, pressure, and vol- 
ume essential experimental variables to characterize the 
relaxation properties. One successful approach to at least 
categorize dynamic properties of supercooled liquids and 
polymers is by expressing them as a function of the ratio 
of mass density, p, to temperature, T, with the former 
raised to a material specific constant 7, viz. 



X = 5(p'^/T) 



(1) 



where x is the dynamic quantity under consideration, 
such as the structural relaxation time r, the viscosity 77, 
or the diffusion coefficient D, and 3 is a function. This 
scaling was first applied to a Lennard-Jones (LJ) fluid, 
with 7 = 4 yielding approximate master curves of the 
reduced "excess" viscosity for different thermodynamic 
conditions [if. More recently Eq. ^ has been shown to 
superpose relaxation times measured by neutron scatter- 
ing y , light scattering [3j , viscosity [4| , and dielectric 
spectroscopy jH, H, 0, HTs] for a broad range of mate- 
rials, including polymer blends and ionic liquids. The 
scaling exponent 7, which varies in the range from 0.13 
to 8.5 PJ3], is a measure of the contribution of density (or 
volume) to the dynamics, relative to that due to tem- 
perature. The only breakdown of the scaling is observed 
for hydrogen-bonded liquids, in which the concentration 
of H-bonds changes with T and P, causing t to deviate 
from Eq. ^ [A]. 

The function 3 in Eq. ([1]) is unknown a priori. Its 
form can be derived from entropy models for the glass 
transitiomleading to an exponential dependence of log r 
on p^/T [m, Another interpretation of the scaling is 
that the supercooled dynamics is governed by activated 



processes with an effective activation energy E{p, T) p^ . 
in which the p-dependence of E{p, T) can be factored 
and expressed in terms of a power law of p. The power 
law scaling arose from the idea that the intermolecular 
potential for liquids can be approximated as a repulsive 
inverse power law (IPL), with the weaker attractive forces 
treated as a spatially- uniform background term 3, 15, 



u(r) 



const 



(2) 



where r is the intermolecular distance. In the case of 
an IPL, in fact, all reduced dynamical quantities [l7j 
can be cast in the form of Eq. ([1]) with 7 = m/3, i.e., 
the thermodynamic scaling is strictly obeyed. For in- 
stance, this applies to the reduced diffusion coefficient 
D* ~ {p^/^T-^/^)D ~ (T3/^-i/2)d. a similar reduction 
of D by macroscopic variables [p and T) has also been 
employed in entropy scaling laws of diffusion ^18'|. 

The IPL approximation emphasizes the dominant role 
of the short-range repulsive interactions for local proper- 
ties such as structural relaxation. Various groups have 
explored through numerical simulations the relationship 
of the steepness of the repulsi ve p otential to properties 
such as the equation of state [3i 0) 21 1, longitudinal 
wave transmission [22J, vibrational spectrum [23i] . liq- 
uid [13] and gaseous|25| transport, the correlation be- 
tween fluctuations of energy and pressure [26| , and the 
fragility [27, 28]. Recently two simulations have appeared 
in which Eq. ([T|) was used to superpose dynamical data 
for polymer chains described using an LJ m-6 potential 
with m = 12 and an added term for the intrachain in- 
teractions. The results appear contradictory: Tsolou et 
al. [i^l obtained a scaling exponent of 7 = 2.8 for the seg- 
mental relaxation times of simulated 1,4-polybutadiene, 
while Budzien et al. ^ superposed diffusion coefficients 
for prototypical polymer chains using 7 = 6 when attrac- 
tion were included in the simulation and 7=12 when they 
were omitted. Thus, the scaling exponent 7 is either less 
than [2§] or greater than [sO] m/3. 
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To clarify this situation and to establish a connec- 
tion between the thermodynamic scaling and the repul- 
sive part of the intermolecular potential, we carried out 
molecular dynamics simulations for supercooled LJ m-n 
liquids, in which the repulsive exponent m was systemat- 
ically varied. Our models are binary mixtures composed 
of = 500 particles enclosed in a cubic box with peri- 
odic boundary conditions and interacting with a LJ m-n 
potential 



Ua(}{r) = 4:ea0 [{(Jap/r)"' - (o-Q/s/r)"] 



(3) 



where a,/3 = 1,2 are indexes of species. We fixed the 
attractive exponent n = 6, as in the standard LJ poten- 
tial, and varied m = 8, 12, 24, 36. The potential Uap{r) 
was smoothed at Tc = 2.bcTaj3 using the cutoff scheme of 
Stoddard and Ford 'si'l. Reduced LJ units are used, as- 
suming (Til, eii and \/mi<7\i/ tu as units of distance, 
energy and time respectively. The mixture on which we 
focus is an additive, equimolar mixture with size ratio 
^ = C22/o'ii = 0.64, equal masses mi = m2 = 1.0 and a 
unique energy scale e^^ = 1.0. The choice m = 12 cor- 
responds to the AMLJ-0.64 mixture studied in H, S^l- 
The samples were quenched isobarically at different pres- 
sures P = 5, 10, 20 by coupling the system to Berendsen 
thermostat and barostat during equilibration (see (s^ 
for details), and performing the production runs in the 
NVE ensemble using the Velocity- Verlet algorithm. The 
timestep 5t was varied according to the repulsive expo- 
nent, ranging from 0.001 (m = 36) to 0.004 (to — 8) at 
high T, and from 0.003 (to = 36) to 0.008 (m = 8) at 
low T . The equilibration criteria were similar to the ones 
used in [s^ . 

The effectiveness of the thermodynamic scaling for LJ 
TO— 6 systems is demonstrated in Fig.[T]for different values 
of the repulsive exponent to. For each to, reduced diffu- 
sion coefficients D* = {p^l^T-^l'^)D were gathered along 
different isobaric paths (P = 5, 10, 20) and the material 
specific scaling exponent 7 was obtained by maximizing 
the overlap between different sets of data, plotted as a 
function of /T. Repeating the analysis for D, instead 
of P**, yields very similar values of 7, but the quality of 
the scaling for D* is slightly superior. Our data span 
roughly 5 decades of variation of P, over about two of 
which the temperature is lower than the so-called onset 
temperature Tq 



34| 



where non-exponential relaxation 
typical of the supercooled regime first becomes appar- 
ent upon cooling the liquid. Analyzing the variation of 
the scaling exponent in our models, we find that 7 in- 
creases with increasing to, but its actual value is system- 
atically larger than to/3. For instance, in the case to = 12 
we obtain 7 = 5.0, a value which we also found to pro- 
vide scaling of D* for other supercooled Lennard-Jones 
(to= 12) mixtures, such as the AMLJ-0.76 mixture intro- 
duced in [sl] and the mixture of Kob and Andersen . 

The origin of the discrepancy between 7 and to/3 lies 
in the fact that the asymptotic region of small interpar- 
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FIG. 1: (color online). Reduced diffusion coefficients D* as a 
function of p^/T for different values of the repulsive exponent 
m at different pressures: P = 5 (squares), P=10 (circles), and 
P = 20 (triangles). From top to bottom: m = 36 (7 = 13.4), 
m = 24 (7 = 9.1), m = 12 (7 = 6.6), and m = 8 (7 = 3.5). The 
estimated uncertainty on 7 is ±0.1 (±0.2 for m = 36). 



tide distances, in which uir) ~ r"™, is not dynamically 
accessible in normal simulation conditions. The presence 
of the fixed attractive term in the potential (Eq. ([3])) 
gives rise to an effective IPL which is steeper in the re- 
gion of r close to the minimum than in the r ^ limit. 
This effect is illustrated in Fig. [2] for the case to = 24. 
The lower panel of Fig. [5] shows a fit of the pair poten- 
tial uii(r) to an IPL (Eq. ^) performed in the range 
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m=24 
P=10, T=0.84 
P=10, T=1.00 
P=10, T=1.50 




TABLE I: Parameters of IPL approximations for Uap(r). The 
effective exponent m is obtained from fitting to Eq. (j4|, 
wfiereas e, k, and x are tlie optimal values for Eq. ((5]). 



FIG. 2: (color online). Top panel: radial distribution func- 
tions between large particles gii{r) at P = 10 for T < To'- 
T = 1.20 (dotted), T=1.00 (dashed), and T = 0.84 (solid). 
Middle panel: gii{r) at the lowest equilibrated T: T = 0.75 
at P = 5 (dotted), T = 0.84 at P=10 (dashed), and T=1.05 
at P = 20 (solid). Bottom panel: pair potential wii(r) (solid) 
and fitted IPL (dotted) in the range [0.95 : 1.01]. The latter 
range is indicated by vertical dotted lines in all panels. 



[ro ■ ri], with ro = 0.95 and ri = 1.01. The value rn — 27.5 
obtained through this procedure is indeed larger than 
m = 24 and is in very good agreement with the value ex- 
pected from the dynamical scaling (37 = 27.3±0.03). The 
range [vq : ri] corresponds to typical distances of closest 
approach between particles probed within our simulation 
conditions, as it can be seen by inspection of the radial 
distribution functions gn{r) (see upper panels of Fig. [2]). 
Extending the range for the fit up to ri = 1.06, which 
is close to the average position of the first peak in the 
(7ii(r), yields a larger value to = 28.8, revealing how 7 is 
dictated by the portion of r around the distance of closest 
approach in the supercooled regime. 

To proceed in a more systematic way, we considered 
all a — (3 pairs (1-1, 1-2 and 2-2) in the potential Uai3{r) 
and performed a simultaneous fit to the following IPL 



J-ai3{r) = e((TQ/3/r)™ + k. 



(4) 
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The range for fitting was defined by two conventional 
distances determined from the radial distribution func- 



tions gapif)'- the distance of closest approach between 
particles, rg, (i.e., the value of r for which the ga/sir) 
first becomes non-zero) and the position correspond- 
ing to half of the height of the first peak, ri, (i.e., 
gapifi) — gapirm)/'^ where Vm is the position of the 
first peak and ro < ri < rm) ■ These quantities depend on 
the thermodynamic state under consideration, but their 
variation with P and T is mild within our simulation con- 
ditions ^38*1. Our interest being the supercooled regime, 
we simply consider the interval [ro : ri] obtained from 
the low-T behavior of the ga/sir)- For each a— /3 pair we 
used the corresponding range [rp : ri] for fitting. In gen- 
eral, the fitted values of m are in good agreement with 
37 (see Table |T| for all values of m. Thus, the scaling 
exponent can be reasonably accounted for in terms of an 
IPL approximation of the pair potential, provided that a 
sensible choice of the relevant range of distances is made. 

The above procedure suggests that a model of soft- 
spheres (SS) with TO = 37 should provide a good refer- 
ence system for the LJ m-6 mixtures. To this aim, we 
approximate Eq. ([3]) with 



t'a,M=r^""/{'^"^^ (5) 
I Ua[3{r) r > xaai3 



where to, e, and k are expressed in terms of x by requir- 
ing continuity of 0th, 1th, and 2th derivatives of Vap (r) 
at r = xaa/3 ■ The value of x is then fixed by requiring 
that 37 = mix) = (toV^"+^ - n^/x"+^)/{m/x"'+'^ - 
n/x"^^). The parameters defining the reference SS mod- 
els for all values of to are reported in TableHl We checked 
that the distance xaa0 always lies in the range [ro : ri] 
defined above. Diffusivity data for the LJ 12-6 mixture 
are compared in Fig. [3] to those of the corresponding ref- 
erence SS mixture along two isochores (/9=1.5, p=1.7), 
which correspond to typical densities attained at low T 
by the LJ system (at constant P). The trend of D{T) for 
the reference system closely follows the one for the full 
LJ system. As expected, the SS mixture has a larger dif- 
fusion coefficient for a given thermodynamic state. The 
contribution to D due to the attractive part of the po- 
tential could also be explicitly included using a WCA-like 
splitting of Vafsir) [1^. For the present purposes, how- 
ever, it is more useful to note that a simple rescaling of 
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FIG. 3: (color online). Arrhenius plot of difltusion coeffi- 
cient D for the LJ 12-6 mixture and the reference SS mix- 
ture (m — 15.0, e — 1.74) along two isochores: p — 1.5 and 
p = 1.7. Inset: reduced diffusion coefficient D* as a function 
of p"^^^ /T. For the SS mixture a reoptimized energy scale 
e — 1.13e was used. 



e (increased by around 10%) yields an excellent super- 
position of D* for all sets of data (see inset of Fig. [3]). 
Thus, at least to a first approximation, the contribution 
of the attractive part of the potential to the dynamics 
alters the shape of the function 3 without affecting the 
scaling exponent 7. 

To summarize, the thermodynamic scaling of the dif- 
fusion coefhcient in supercooled LJ m-Q mixtures reflects 
the importance of the repulsive part of the pair potential 
in determining the dynamical properties of these systems. 
The scaling exponent 7 is larger than m/3 for LJ m-6 
liquids, a fact which can be rationalized by approximat- 
ing the repulsive part of the potential with an IPL having 
exponent m w 87. Generalizing such arguments to more 
realistic models of glass-formers 2^ 3^ and establishing 



connections with other scaling procedures for D* isl. 37 1 
are open challenges for future investigations. 
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